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The s = 1 spinor Bose condensate at zero temperature supports ferromagnetic and polar phases 
that combine magnetic and superfluid ordering. We investigate the formation of magnetic domains 
at finite temperature and magnetic field in two dimensions in an optical trap. We study the general 
ground state phase diagram of a spin-1 system and focus on a phase that has a magnetic Ising 
order parameter and numerically determine the nature of the finite temperature superfluid and 
magnetic phase transitions. We then study three different dynamical models: model A, which 
has no conserved quantities, model F, which has a conserved second sound mode and the Gross- 
Pitaevskii (GP) equation which has a conserved density and magnetization. We find the dynamic 
critical exponent to be the same for models A and F (2 = 2) but different for GP (z « 3). Externally 
imposed magnetization conservation in models A and F yields the value z ~ 3, which demonstrates 
that the only conserved density relevant to domain formation is the magnetization density. 

PACS numbers: 03.75.Mn, 03.75.Kk, 64.60.Ht, 75.40.Gb 



INTRODUCTION 



The field of cold atomic gases has witnessed an explo- 
sion of experimental and theoretical research in the last 
ten years. The study of these systems has combined ideas 
from various disciplines of physics such as atomic physics, 
condensed-matter physics, optics etc. Cold atomic sys- 
tems have provided a testing ground for some of the most 
fundamental principles of collective quantum behavior 
like Bose-Einstein Condensation. Of particular interest 
is the study of spinor condensates, which are condensates 
of atoms with non-zero spin and have been the focus of 
intense experimental^i^i^i^ and theoretical^i^i^ii^ studies in 
recent years. The spin degree of freedom opens up the 
possibility of interesting collective magnetic behavior in 
these systems in addition to the phenomenon of Bose- 
Einstein condensation. It has already been demonstrated 
that the presence of spin greatly modifies the nature of 
the condensate and superfluid transition in spinor con- 
densates compared to those without spiii^i^. 

Spinor condensates have over the last few years been 
realized in both magnetic and optical traps. The lat- 
ter are more interesting from the point of view of spin 
ordering, since the spin degree of freedom is not frozen 
out. The most widely studied atomic systems are those 
of the spin-1 alkali atoms ^^Na and *^Rb. These sys- 
tems differ from each other in the nature of the effective 
two-body interaction, which is antiferromagnetic in the 
former and ferromagnetic in the latter. The condensates 
with antiferromagnetic interactions are also called polar. 
Recent advances have made it possible to image ferro- 
magnetic domains in optical traps of *^Rb, making it 
possible to study the interesting physics of domain forma- 



tion in them''. This technique requires the application of 
a magnetic field, an additional tunable parameter which 
makes the phase diagram of these systems interesting. 
Moreover, these atoms have also been trapped in two di- 
mensional geometries, where the physics of collective be- 
havior is often more exotic than in higher dimensions^i^. 
The importance of this experiment for basic condensed 
matter physics is twofold: it probes both our understand- 
ing of phase-ordering kinetics at finite temperature (when 
observed at the longest times) and, as the temperature is 
lowered or the observation time is shortened, our under- 
standing of dynamics across quantum phase transitions. 

In this paper we will investigate magnetic domain for- 
mation in spin-1 systems at finite temperature and mag- 
netic field. The main purpose of this study is to compare 
and contrast various plausible dynamical models with re- 
spect to coarsening of a magnetic order parameter. The 
quantity of primary interest, will be the dynamic critical 
exponent z which determines the rate of domain forma- 
tion at large times: the domain size L grows with time 
as L . We will examine the general phase diagram 
of spin-1 condensates in the presence of a magnetic field 
in an optical trap and comment on the broken symme- 
tries of the various ordered phases. We will then choose 
the phase that is most convenient to a study of magnetic 
domain formation and elucidate the similarities and dif- 
ferences between dynamic models, highlighting the im- 
portance of different conservation laws in the dynamics. 
We will compare our results with existing ones wherever 
possible. 

A natural question is how the stochastic time- 
dependent Ginzburg-Landau (TDGL) approach in this 
paper is related to previous studies using deterministic 
equations of motion, such as the Gross-Pitaevskii equa- 



tion for the condensate, plus quantum kinetic theory for 
excited statea^ '^^i"'^"'^ . The answer is that the correct de- 
scription depends on experimental parameters such as 
the time scale of observation and the normal-state pop- 
ulation. The time scale at which stochastic processes 
resulting from interaction with the normal cloud become 
important can be increased by decreasing the tempera- 
ture of the system. The initial instability in a finite trap 
is likely to be described correctly by the deterministic 
theories in the literature; coupling to the many degrees 
of freedom in the normal cloud is irrelevant for the imme- 
diate dynamics of the condensate. However, the longer 
times accessed in current and future experiments are ex- 
pected to be described by the theory developed here. In 
other words, the universal dynamical properties in the 
sense of critical phenomena are described by the theo- 
ries presented here at any finite temperature, as long as 
the system is observed for a sufficiently long time. We 
believe that current experiments may already be in the 
regime where the theory presented here is valid. How- 
ever, even if they are not, increases in observation time 
will soon enable a precise comparison between theory and 
experiment. 

Our main results on phase ordering of spinor conden- 
sates are contained in sections VH and VHI. We argue 
in the final discussion that one specific dynamical model 
("model F" dynamics, in the notation of the review paper 
of Hohenberg and HalperirtU') is expected to describe the 
long-time dynamics of spinor condensates. This dynami- 
cal model is a more complicated version of the model used 
in earlier studies of superfluid o^^'^'^i^^'^^ , and reproduces 
the known propagating modes of the spinor condensate 
at zero temperature. All parameters in the dynamical 
model can be determined from measurements of the con- 
densate, as explained in the appendix. 



II. THE MAGNETIC PHASE DIAGRAM OF 
SPIN-1 BOSONS IN AN OPTICAL TRAP 



Spin- 1 condensates are theoretically more complex 
than those with zero spin^-''^^ in that the condensate order 
parameter is a three component complex vector 



* = 




(1) 



with tpa being the order parameter in the spin state of 
eigenvalue a along some arbitrarily chosen direction. If 
one assumes that the condensate state is a single particle 
zero momentum state, the total energy for a given density 
of atoms in an optical trap with a magnetic field B in the 
z direction can be written as 



Here S = SxX: + Syij + SzZ, where 
Sx — 




(3) 



are the generators of SU{2) in the spin-1 representa- 
tion and (A) = 'i'^'A'ii. ci is the spin-spin interaction 
which can be antifcrromagnetic [ci > 0) or ferromag- 
netic (c2 < 0). g2 oc and the second term is just the 
quadratic Zeeman term. The absence of a linear term is 
due to the fact that the time for the relaxation of mag- 
netization in optical traps is less than the lifetime of the 
condensate itself. The ground state manifolds can be 
obtained by minimizing the free energy with respect to 
{tpa}- It already been shown that in the absence of 
a magnetic field, the ground state manifolds in the po- 
lar and ferromagnetic cases are isomorphic to the spaces 
'^^^^^'^ and 50(3) respectively- . The phase diagram in 
the presence of a magnetic field is given below. 
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FIG. 1: The ground state phase diagram of a spin-1 conden- 
sate in an optical trap in the presence of a magnetic field 
that couples through a quadratic Zeeman term. The different 
quadrants have different phases with various types of in-plane 
and out-of-plane ordering. This figure has been taken from 
Mukerjee et. a&. 

Fig. [T] has four quadrants labelled by signs of C2 and 
(72- In the polar (c2 > 0) case, the magnetic ordering is 
either in plane or out of plane depending on the sign of 
(72- "Magnetic ordering" here refers to the ordering of the 
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spin-quantization axis (fi). The ground state is always a 
macroscopically occupied single particle state of zero spin 
projection in this case. For 52 > 0, the only symmetry 
that is broken in the ordered state is that of the U{1) 
phase {9) of the condensate. For 172 < 0, however there 
is an additional U{1) due to the in-plane ordering of the 
spin-quantization axis. The phase and spin are coupled 
through a Z2 identification, which denotes symmetry un- 
der 9 ^ 9 + n and n — s- — n. The vortices corresponding 
to the spin and phase are thus coupled and can lead to in- 
teresting finite-temperature physics in two dimensions'^. 

The lower part of the phase diagram corresponds to the 
ferromagnetic case (52 < 0) and will be of primary inter- 
est to us. The lower left quadrant corresponds to the case 
(72 < 0. The ground state now breaks a, U{1) symmetry 
corresponding to the phase and an Ising Z2 symmetry 
corresponding to the spin. Physically this means that in 
the condensate, the bosons are either in a state of spin 
projection 1 or -1. It is this Ising degree of freedom, 
we exploit to study domain formation in two dimensions. 
The reason is that since long range Ising order is possi- 
ble in two dimensions (as opposed to U{1) order), it is 
easier to define and measure the sizes of large magnetic 
domains required to investigate long time behavior. It is 
thus this quadrant that will be the focus of the rest of 
out studies. For the sake of completeness we note that 
the lower right quadrant, which corresponds to the case 
(72 > is divided into two parts by a straight line with 
equation 52 = 2c2. To the left of this line, one has in- 
plane ferromagnetic ordering with the spins pointing in 
some U{1) direction in the plane. The ground state thus 
breaks two U{1) symmetries, one corresponding to the 
phase and the other corresponding to the spin. To the 
right of the line, it is energetically favorable for the sys- 
tem to be in a polar out of plane state. This suggests 
the interesting possibility of a quantum phase transition 
in these systems tuned by the magnetic field. 



III. 2D FINITE TEMPERATURE PHASE 
TRANSITIONS 



Since we are interested in studying finite temperature 
coarsening dynamics in the 2D system with C2 < and 
g2 < 0, it is important for us to locate the position of 
the superfiuid and magnetic transitions. This problem is 
also interesting in its own right since such situations also 
come up in the study of classical frustrated spin systems, 
like the fully frustrated XY antiferromagnet (with tt flux 
per plaquette) on a square lattice or the triangular lattice 
XY antiferromagnet, where the Z2 corresponds to a chi- 
rality. The U{1) and Z2 transitions are in close proximity 
to each other in these cases. The situation is our partic- 
ular case is not very different. We find that the U{1) 
transition is of the Kosterlitz-Thouless (KT) type and 
the Z2 transition of the 2D Ising type. Furthermore, we 
find that for a certain range of parameters, T^^ > 11/(1) 
which is also what is observed in the fully frustrated XY 
model on the square latticoi^ and for others the order 
of the transitions appears to be reversed. This depends 
on the magnitude of the ratio of the parameters §2/02- 
For small values of this ratio, T^^ > T^(^ly There is pre- 
sumably also a point where the two transitions occur at 
exactly the same temperature, where the combined tran- 
sition is in a different universality class from 2D Ising 
and KT. We present here numerical data on just one set 
of parameters where > Tkt and illustrate how the 
two transitions can be accurately determined despite be- 
ing reasonably close to each other in temperature. The 
method used is due to Olsson^^ and we employ a nu- 
merical Monte-Carlo simulations that uses the following 
Ginzburg-Landau free energy functional 



J 



F ^ dr 



C2 



(4) 



r 



with the following set of parameters, {a = 0.5, = 
5.5, CO = 7.0, C2 = -2.4,52 = -1.3}. 

The Kosterlitz-Thouless (KT) transition is detected 
by observing the temperature dependence of the helicity 
modulus Y . The helicity modulus for a discrete system 
of N lattice points is defined as 



r = 



1 < > 



N dS^ 



(5) 



5=0 



where 5 is a fiux twist applied along a particular direc- 
tion. The helicity modulus undergoes a jump of magni- 
tude at the location of the transition. This is shown 
in Fig. [21 where the transition temperature is seen to be 
Tc w 0.44. 

The standard method to determine the location of the 
Z2 transition using fourth order cumulants of the mag- 
netic order parameter fails here for the same reason that 
it does in the case of the fully frustrated XY model, 
which is the proximity to the KT transitionii. The cu- 
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FIG. 2: The helicity modulus as function of F for the param- 
eter set {} for different system sizes. A clear jump is visible 
of 2Ta/-K is visible at Tc ~ 0.44. 



and determine Tc by fitting it to the expected 2D Ising 
form. The magnetization M{r) is given by 

M(r) = |V^+i(r)|2-|V._i(r)|2. (6) 

The correlation length ^{T) can be extracted from the 
magnetic autocorrelation function 

g{r) = (M(r)M(O)) = e-^'^^^\ (7) 
If the transition is 2D Ising like, 

C(T) ^ (8) 

The numerical result is shown is Fig. [31 which shows that 
the correlation length fits the 2D Ising form fairly well. 
The obtained transition temperature is ~ 0.53. A 
more careful finite-size scaling analysis can be done to 
determine the two transition temperatures, but even at 
this level of analysis it is clear that T^.-, > 



IV. DYNAMICAL MODELS 
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FIG. 3: The magnetic correlation length as a function of tem- 
perature fitted to the 2D Ising form. Tc ~ 0.53 as estimated 
this way. 

mulant method assumes that the only relevant length 
scale at the transition is the system size which is not 
true here because of the large correlation length corre- 
sponding to the closely situated KT transition. Thus, a 
more accurate method is to look at the critical exponent 
of the correlation length of the magnetic order parameter 



The study of the formation of domains of the order pa- 
rameter requires careful consideration of the dynamical 
modes of the system. Dynamical models are often con- 
strained by conservation laws that are present as a con- 
sequence of symmetries or otherwise in the system. It is 
well known that the presence of conservation laws usually 
affects the rate of formation of domains, since the phase 
space of states that the system can pass through in the 
approach to the ordered state is constrained by the con- 
servation laws. However, not all conservation laws affect 
domain formation in the same way and some might be 
more important than others. In this section we consider 
some dynamical models appropriate for the description 
of our system and comment on the conservation laws and 
the dynamical modes obtained from them. 

The most commonly used dynamical model to describe 
spinor condensates is the Gross-Pitaevskii (GP) equa- 
tion"'^^. This is the model that has been extensively used 
to study domain formation in these systems. The model 
consists of treating the condensate as a classical field 
at zero temperature whose dynamics are given by the 
Hamilton equations of motion of the appropriate Hamil- 
tonian. In our case, the Hamiltonian is 



H = dr 



(9) 
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with the dynamical equation of motion 



ih 



dt 



SH 



(10) 



It should be noted that the condensate density V'aV'a 
and magnetization are both conserved by this equation. 
However, the GP equation cannot correctly describe the 
approach to equilibrium at finite temperatures, since 
the dynamics is only precessional and not relaxational. 
While this equation might be appropriate for the de- 
scription of the dynamics once the condensate has been 
formed, it is inappropriate for the study of dynamic phe- 
nomena in other cases, for example quenches from high 
temperature where the energy of the condensate is not 
conserved. 

The effect of finite temperature on the dynamics in 
spinor systems has thus far been taken into account 
through phenomenological rate equations, which too do 
not describe the approach to equilibrium. A simple 
model which is more appropriate is the so-called "model 
A" of the Hohenberg and Halperin classification^^. This 
model uses the Ginzburg-Landau free energy Eqn.[4]with 
simple Langevin dynamics. Operationally, this means 
the dynamical equation 



dt 



-ro 



SF 

Ha 



Ca(r,t). 



(11) 



Thermal fluctuations due to finite temperature are con- 
tained in the noise variable Ca(r, t), which has the follow- 
ing autocorrelation function 

{C{r,t)Cbir',t')) ^2Re{To)kBTS{r-r')dit-t')dab (12) 

consistent with the fluctuation-dissipation theorem, that 
drives the system to equilibrium from a non-equilibrium 
state. This model like the GP model only considers the 
condensate as a classical field but unlike the GP model 
does not possess any conservation laws. It is relaxational 
in nature with the rate of relaxation of the order param- 
eter set by Re Fq and can be a complex number. The 
condensate density is no longer conserved and neither is 
the magnetization. The condensate is exchanging parti- 
cles and energy with the "normal fluid" in this model. 
The non-conservation of energy of this model can be rec- 
tifled by implicitly including the normal fluid through 
a conserved "second sound" mode (m), which is a real 
scalar field. Based on previous experience with the su- 
perfluid transition in helium, one expects that a correct 
description of the dynamics near the transition requires 
this additional field and in the notation of Hohenberg 
and Halperin, this model "model F" , with random forces 
and 6 for the fields V'a and m respectively. The free 
energy i^gg with this second-sound mode is 



where F is given by Eqn. [J] The dynamics are given by 

(14) 



at dyj* dm 



Fss = F - 



(13) 



IP = A?rV^-- + 2,oInr^^:--J+.(r,t) 

These equations conserve the second sound density m 
and the noise correlator for 0, consistent with the 
fluctuation-dissipation theorem is 

(r(r, t)T(r', i')) - -2A™V25(r - v')E(t - t!) (15) 

The free energy _Fgg has extra terms compared to F\ a 
term that couples m and the condensate density and 
others that contain the energy of the mode m. The 
dynamical equations also incuding coupling terms as 
a consequence of the non-vanishing Poisson brackets 
{to, V'a] Aside from terms that result from deriva- 

tives of the Ginzburg-Landau free energy, there could 
in principle be additional magnetic terms analogous to 
those in the Heisenberg ferromagnet (e.g., Sx V^S, where 
S is the local spin density). Since this is even in S, it 
will not be obtained as the S derivative of any free en- 
ergy, but will originate in the microscopic Hamiltonian. 
A check that no such additional terms are necessary is 
that the above equations reproduce the previously ob- 
tained modes in the GP equation. Such a calculation 
was carried out by Hohenberg and Halperin for the case 
of Helium and we extend that to the case of spinor con- 
densates in the next section. 

The three different dynamical models, GP, model A 
and model F are the ones we will use to investigate do- 
main formation in spinor condensates at flnite temper- 
ature and magnetic flelds. We will in addition to the 
dynamical equations above also impose the conservation 
of magnetization on models A and F, to investigate the 
effect of that conservation law on the dynamics. As can 
be seen from the above discussion, model F contains 
many more parameters than model A and the GP equa- 
tion. The parameters of this model are related to possible 
experimentally-measurable quantities in the appendix. 



DYNAMICAL MODES IN THE ORDERED 
STATE 



Model F contains in it both the GP equation and model 
A, which can be seen by setting the appropriate param- 
eters in it to zero. However, it is important that model 
F produces all the dynamical modes that the GP equa- 
tion does even when these parameters are not zero, in 
order for this treatment to be valid. There will also be 
additional modes produced (for example in to), that are 
absent in the GP equation. We explicitly demonstrate 
this in this section. 
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We begin by setting the temperature and magnetic 
field to zero, to enable comparison with the GP equa- 
tion. The idea is to check that the introduction of the 
extra parameters of model F does not alter the modes 
that have already been calculated — . It is known that 
there are three linearly dispersing mode in the polar case 
and one gapped, linear and quadratic mode each in the 
ferromagnetic case. The dynamical equations describing 
the modes (either propagating or diffusion) are 



dt 

dm 

'dt 



-ro 



5Fss 
Hi 
dFss 
dm 



SFss 



(16) 



2.oIm(C^) 



For brevity of notation, we also set T — 
aoT^^ and exphcitly write Fo = Fi + iFs 
Fi and F2 are real. 



A. The polar case 



-aV^, /i = 
where both 



We assume that 



where 





N = ^\ 1 




(17) 
(18) 



is the value of the order parameter in the ordered polar 
state and $ is a perturbation on it. The number of parti- 
cles within the condensate, is related to the value fi and 
Co by minimizing the free energy Fgs 



(19) 



Notice that in order to get these equations, we need to 
take exactly the value in (fT9|) ) 

The equations for (jji and 0_i have the same form as 
for the GP equation^, so the spin wave modes are the 
same. Both = 0+1 + i<j)-i and 7\f_ = (p+i — 
disperse linearly with k, with velocity cm — ^2\/noC2. 

The step by step solution for the coupled equation be- 
tween (jjQ, 4>Q and m is tedious, so we only write down 
the result here. Basically, the second sound mode and 
density fluctuation Shq = 0q -t- (/iq couple and form two 
modes with linear dispersion relations, the velocity is 



C 



2aT2ConQ 



(21) 



Notice that the first term in the square root in the 
above equation is the square of the second sound veloc- 
itjii^, with r2 = and ignoring the propagating mode of 
■0. The second term in the square root is the one that 
appears as the density fluctuation mode in Ref. 5, where 
the the second sound mode was ignored. Here we see that 
if we take into account both densities, the second sound 
mode and the density fluctuation mode couple into a new 
mode with velocity Cg. 



2. Polar state with Fi = 0, Ao / 0, 7 = 



Here the dissipation Aq term is added back into the 
equations. We will not consider the case with finite Fi, 
since we assume that within the condensate, the dissipa- 
tion of the modes is very small. 



The coupled equations between m and 00 now become 



1. Polar state with Fi = 0, Ao = , 7 = 



Let us first ignore the coupling jmip^ipa, as well as all 
the dissipation terms in the equations, like the term with 
coefficient Fi and Ao. We will put them back in later. 

In this case, all the modes are propagating, since there 
is no dissipation. After expanding the equations around 
iV, the linearized equations we have are 



-«.9o\M)V^(0o - (f>o) 
-iF2(T0i + noC2(0i + 01i)) 
-iF2(T0_i+noC2(0^ +0-i)) (20) 



= iF2(r0o + 2cono(0o + 0o)) + *"^-^5o 

dm 
'dt 

dt 
¥-1 
dt 



^ - ^F2(^0o + 2cono(0o + 0S))+^m^5o(22) 



dm — 2/ 

-Q^ = -«.goV"oV ( 



Ao„2 



(23) 



The detailed solution is again tedious and we will solve 
the equation based on following approximation that the 
higher order terms of spatial derivatives are small, since 
we are only interested in the limit k going to zero. Under 
this approximation the dispersion relation is 



LJ — Cgk + i-^k^ 



(24) 



The mode gets a propagating part, which is linear in 
fc, and a damping part, which is proportional to k^. Cg is 
given by (pij) . 
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Polar state with Fi = 0, Ao 7^ 0, 7 7^ 



Turning on 7, changes only two terms in the equations. 
First, 



dt 



-iT2T(/)_i + 202^00-1 



(31) 



The dispersion relation for the second sound mode is 



= ir2{T(l)o+2cono{(l)o+(l)Q))+igojno{(j)Q+<j>^)+im—^gQ. 



We can redefine 



Co = Co + — 
i 2 



(25) 



(26) 



and make the equation look exactly like ([20]), except for 
replacing co by Cg . 7 also modifies (|23p , by adding a term 
to the right hand side 



— — = ir2(T(j)o + 2cano{(j)o + 4>o)) + im 



dt 
dm 

'dt 



C 



-go 



-^V^v2(0o + 0S) 



(27) 



Solving this modified equation, we see that the term pro- 
portional to 7 only contributes higher order momentum 
terms. So, up to linear order in fc, the dispersion rela- 
tion is not changed. Therefore all the results here are the 
same as case 2, if we replace co by Cg. 



B. The ferromagnetic case 



(32) 



C 



2ar|(cg + C2)ng, 



where we have kept terms only to lowest order in the 
momentum in every step of the calcultion. For the den- 
sity fluctuation mode Sn = ^/no((j)i + (f)i'), the dispersion 
relation is 

u! = Cgk (33) 

The spin wave mode 5M_ — y^rio(j)Q has the same dis- 
persion relation as the one obtained in the GP case^, 
Lu = ak^. Again, turning on the interaction 7 does not 
change the result much. It causes a redefinition of cg in 
the propagating part of the mode, and only contributes 
higher order momentum terms in the diffusion or damp- 
ing part of the mode, which are not important when the 
momentum is small. 

Thus, we see that in both the polar and ferromagnetic 
cases, the dynamic modes obtained in the presence of the 
extra parameters of model F are consistent with those ob- 
tained from the GP equations. The nature of the density 
mode changes because of coupling with the second sound 
mode but the spin wave modes remain unaffected. 



The solution of the ferromagnetic case is very similar 
to the polar case. The general formalism and effective 
action Eqn. [13] still apply. The difference is in how we 
linearize the equations. In the ferromagnetic case, we 
should linearize the equations around the state 



^ = N + 



with 



N 




(28) 



(29) 



Here the density of the condensate is not only related to 
the coefficient cg, but also to the coefficient C2. 



" Cg -I- C2 

We now obtain the following linearized equations 



(30) 



= iT2{T4)o + 2{co + C2)no{(f)i + tfiD) + im^^go 



dm 
'dt 



VI. DOMAIN FORMATION 



A typical experiment or numerical simulation of 
coarsening involves starting the system off at a high- 
temperature (usually disordered) state and rapidly 
quenching it to a temperature below the ordering tran- 
sition to observe the growth of domains of the ordered 
state. At the heart of the theoretical analysis of this 
process is the scaling hypothesis2^. The equal-time cor- 
relation function of the order parameter m(r, t) is defined 
as 



C(v,t) = (m(x + r,t)rn(r, i)). 
The scaling hypothesis states that 

C(r,i) = / 



L{t) 



(34) 



(35) 



L{t) is a characteristic length scale, the domain size. Fur- 
ther, at long times L cx i^^^, where 2, the dynamical 
critical exponent. The dynamical critical exponent is 
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dependent on the model used to describe the ordering 
dynamics of the system, the symmetry of the order pa- 
rameter and the nature of defects present in the initial 
state. For instance, it is known that z — 2 for the Ising 
model with dynamics which do not conserve the total 
magnetization, after a high temperature quench. For an 
XY model on the other hand z = 2 but with a loga- 
rithmic correction as a function of time^^,. The difference 
from the Ising case can be attributed to the different bro- 
ken symmetry and consequently the topological defects 
present in the high temperature state. The Ising model 
with a conserved order parameter on the other hand pro- 
duces a different dynamic critical exponent z = S^i^i. 
The growth of domains is slower in this case compared 
to the case with no magnetization conservation since the 
conservation law places constraints on the phase space 
available during domain growth. 



VII. DETAILS OF THE NUMERICAL 
SIMULATION 



In this work we study the dynamics of domain growth 
in 2D after a magnetic field quench and not a temper- 
ature quench. The motivation is a similar approach 
adopted in recent experiments in optical traps. To be 
specific, we study a ferromagnetic condensate in two di- 
mensions whose initial state is a polar out-of-plane state 
in the fourth quadrant of Fig. [1] and quench it to a value of 
the field, where the ordered state is a ferromagnetic out- 
of-plane state (in quadrant 3 of Fig. [T]) . Operationally, 
we sweep the parameter 172 from a large positive value 
to a negative value. For models A and F, this is done 
at finite temperature and the order parameter eventually 
relaxes to a uniform value consistent with the ferromag- 
netic out-of-plane state. The Gross-Pitaevskii equation 
on the other hand does not cause the system to relax 
but rather to oscillate between different concentrations 
of the three spinor components. Further, it does not al- 
low an initial state which is completely polar out-of-plane 
to form magnetic domains of the -1-1 and -1 components 
at any value of the time. In this case, we start with an 
initial state, which has 90% of the atoms in the (polar 
out-of-plane) state and the other 10%, divided equally 
among the -1-1 and -1 states. The phase of each spinor 
component in the initial state is chosen to be a random 
number between and 2tt allowing for spatial inhomo- 
geneity which leads to domain formation. 

The equations of motion corresponding to each model 
are integrated numerically using a first order Euler 
method with the noise functions drawn from a Gaussian 
distribution. The size of the numerical grid ranges from 
50 X 50 to 200 X 200. The time step is adjusted depend- 
ing on the values of the other parameters and varied to 
check for consistency. The number of parameters is large 



(especially for model F) and we present results only for 
a fixed set of parameters. However, we have explored 
other parts of the parameter space consistent with fer- 
romagnetic out-of-plane order and not found any qual- 
itative and wherever appropriate (like for the value of 
z) quantitative difference in the results. The set of pa- 
rameters for which we report results are those in section 
III with the additional model F parameters, {Re(ro) = 
1.30,Im(ro) = 0.26,30 = 0.35, Ag" = 0.84, 70 = 1.5}, 
wherever applicable. 

The domain size L[t) is measured using the relation 

where So{t) and S2{t) are respectively the zeroth and 
second moment of the structure function 

S'(k,<) = y"dr(M(r,^)A/(0,^))e*•^ (37) 

which is the Fourier transform of the order parameter 
correlation function. The domain size is also calculated 
by measuring the size of domain boundaries directly in 
the simulation grid. This method serves as a consistency 
check on the first method. It should be mentioned though 
that the second method is useful and consistent with the 
first one only when there are very few small bubbles of 
one value of the order parameter inside large islands of 
the other value. This method essentially ignores these 
bubbles by looking for large closed domain walls and 
works best when the domains are large in size. 



VIII. RESULTS 



A. Finite temperature without conservation of 
magnetization 



We present results for the domain size as a function of 
time for models A and F in Fig ID The results presented 
are for the set of parameters mentioned in the preceding 
section, with a magnetic field quench and have been ob- 
tained on a grid of size 200 x 200. It can be seen that 
domain formation is faster for model A than for model F, 
which can be attributed to the presence of the extra con- 
servation law. This certainly appears to be the case over 
the range of parameters that we have explored, but may 
not be the case elsewhere in parameter space. Whether 
or not this is a universal feature requires more careful 
analysis. The curve for Lit) as a function of t for model 
A dynamics seems to yield a 2; = 2 ± 0.15 over the entire 
range of values of time we have presented. Further, the 
value of z obtained at different values of time seems to 
be fairly constant. This is the value of z, one would ex- 
pect for a high temperature quench in a pure Ising model. 
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to first order in the domain size 



* 



* 



■ model A 
* model F 



z(t) zit ~ oo) 



1 - 



m 



(38) 



This suggests that z{t) approaches its infinite time value 
from above, which appears to be the case here as well as 
can be seen from Fig. [5] , which is a plot of l/z{t) vs. 
1/L(t). However, we emphasize that the above analy- 
sis is strictly applicable only to the case where the or- 
der parameter is conserved, which is not the case here. 
The quantity that is conserved is the second sound mode. 
Nevertheless, it is possible that the drift can be explained 
by some mechanism similar to the above. 

To conclude this part, we remark that both models A 
and F without any explicit magnetization conservation 
both yield the same dynamic critical exponent z = 2 for 
coarsening with a magnetic field quench. 



FIG. 4: L{t) as a function of logjQ t for model F with the 
parameter set TZ, with no magnetization conservation . 



B. Finite temperature with conserved 
magnetization 



0.46 
0.1 



0.2 



0.3 
1/L(t) 



0.4 



FIG. 5: l/z{t) as a function of 1/L{t) for model F and no 
magnetization conservation with the parameter set TZ demon- 
strating the drift of z{t) as a function of t and thus increasing 
L{t) towards the value z = 2. 
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FIG. 6: L{t) as a function of log^Q t for models A and F and 
conserved magnetization density with the parameter set TZ. 



Model F also yields z ^ 2. Unlike in model A dynam- 
ics, there is a small drift in the value of z obtained at 
different values of t. A similar drift (of a larger magni- 
tude) has been seen in the case of the Ising model with 
dynamics that conserve magnetization and it has been 
argued by Huse^** that this is due to excess transport in 
domain interfaces. It then follows that the effective dy- 
namic critical exponent z{t) drifts in the following way 



We now present results for models A and F with con- 
served magnetization. The magnetization conservation 
is implemented in terms of a local continuity equation 
in the magnetization density and a magnetization cur- 
rent. We illustrate how we do this for model A and 
the implementation for model F proceeds along simi- 
lar lines. We first note that the magnetization density 
M — |V'+i(r, i)p — |'0-i(r,i)P only involves the ampli- 
tudes of the componenets of the condensate order param- 
eter. We first write down model A dynamics in terms of 
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FIG. 7: Ijzit) as a function of for models A and F 

and conserved magnetization density with the parameter set 
IZ demonstrating the drift of as a function of t and thus 
increasing L(i) towards the value z = 3. 



separate dynamical equations for the phase and ampli- 
tude of each component of the condensate order param- 
eter. These turn out to be 



3Q 1 1 

(40) 



where the noise correlators are 

(/ia(r,i)/xb(r',i')) = (y^(r,t)v^,(v\t')) (41) 
= Re(ro)fcBT5(r-r')5(t-0'5afc. 

Note that every quantity in the above equations is now 
real. If we were interested in conserving the density |'0aP 
of each component individually, we would modify Eqn. l39l 
to 



dt 



im(ro)v^ ( i; 



+Aia(r,0, (42) 
with the correlator for /i now given by 

(Ma(r,t)Ma(r',i')> =4Re(ro)V2[|^,(r,t)p5(r-r')]5(t-t')- 

(43) 

This ensures there is a conservation equation of the sort 



dt 



-V.Ja(r, t) 



(44) 



2|^a 



1 5F 

Mro)72r + A^'«(i"'0' 



(39) 



for each component. We are however not interested in 
conserving the density of each component, but only the 
combination M = — IV'-iP- To this end, proceed- 

ing as above, we obtain the following set of equations. 



dM ^ , , 9 f, , , SF , , , \ , , ^ f 5F SF 



dN / SF SF \ ( SF SF \ 

Here N = Iv^+ip + and the noise correlators are given by 

{fiM (r, t)fiM (r', t')) = 4Re(ro) V2 [{ | V^+i (r, t)\' + \yj+,{r, t)\^}S{r - r')]S{t - t') 

{fiNir, i)MJv(r', t')) - 4Re(ro){|^+i(r, t)\^ + |V^+i(r, t)\^}S{r - r')S{t - t') 

(Aio(r,i)Mo(r',0> = B.e{T^)kBTS{v - v')S{t - t') 

{va{Y,t)yb{v\t')) = Re{Vo)kBTS{v-Y')S{t-t')Sab (49) 



r 



The noise functions hn, /^a/, Mo and Va are mutually un- 



correlated. The above equations for M and N can be 
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used to generate equations for and which 

is the way the numerical calculation is performed. Note 
that the full set of dynamical equations written above 
has no conservation law except the one for M. We now 



use the same procedure to impose magnetization conser- 
vation on model F. The dynamical equations in this case 
are 



dM 
~dt 
dN 
'dt 

dt 

dm 
'dt 



-Re(ro)V2 j^lV'+i 
-Re(ro) f iV^+i 



(5|-0+i 



1 



■ - 1^-1 



+ Im(ro)V 
Im(ro) 



'SFss 


SFss 




S9-1 


SFss ^ 




69-1 _ 





1r ^ 



SF 

Im(ro)^+Aio(r,t) 



1 p .J. JFss ^ U .SFss 
A^V^^ + onlml'l^J^^ 



dF 

90\i^a\-Q^+Va{v,t) 



6m 



+ HM{r,t) 



/lAr(r, t) 



(50) 
(51) 
(52) 
(53) 
(54) 



The noise correlators are the same as for model A with 
magnetization conservation with the additional correla- 
tor 

(r(r, t)r(r', t')) - -2X^WH{r - r')6{t - t'), (55) 

as in the case of model F without magnetization con- 
servation. It should be noted that the coefficient go ap- 
pears only in the dynamical equation for the phases of 
the different components of the condensate thus making 
its identification as a precessional term obvious. Further, 
the above dynamical equations conserve both the mag- 
netization M and the second sound mode m and thus 
represent perhaps the most realistic dynamical model for 
a BEC at finite temperature and field; one where the con- 
densate can exchange charge and energy with the "nor- 
mal cloud" but not magnetization. 

Once again, the results presented are for a 200 x 200 
simulation grid. We have checked that the additional 
magnetization conservation law does not affect the static 
properties of the model. As in the previous case, we 
again see that domain formation is faster for model A 
and than model F. This time, however, the dynamic crit- 
ical exponent obtained is not equal to 2. As can be seen 
from Fig. [3 which is a plot of l/z{t) vs. 1/L{t) for both 
models, there is a significant drift of z(t) as a function of 
time. Once again the direction of the drift is consistent 
with Eqn. I38l and in this case, the analysis mentioned in 
the previous subsection is directly applicable, since it is 
the order parameter (the magnetization) that is directly 
conserved. However, the noise levels of the simulations 
do not permit a fit to Eq. [38l It should be noted though 
that to the extent observable in the numerical simula- 
tion, there is a drift in the direction of z = 3 in the data. 
This is what is observed in a pure Ising model with con- 



served magnetization in a high temperature quench. One 
interesting observation is that the value of l/z{t) seems 
to deviate more strongly from Eq. [35] at large times for 
model F than model A. This deviation has also been ob- 
served for the simple Ising model with conserved magne- 
tization, where it was attributed to finite-size effects and 
correlated noise in the simulations. That could well be 
the case here as well, although it is not clear why these 
effects should be more pronounced in one model than in 
the other. It should be noted that in the simulations on 
the Ising model^i, the values of z{t) observed for com- 
parable simulation times are roughly close to what we 
observe. 

The conclusion of this part is that models A and F with 
explicit magnetization conservation yield a dynamic crit- 
ical exponent z « 3, different from the exponent obtained 
without magnetization conservation. Thus, the models A 
and F seem to be identical as far as long-time coarsen- 
ing behavior of the magnetization is concerned and the 
behavior is truly determined by whether or not the mag- 
netization is conserved, which it is not explicitly in either 
model. The difference between these two models will be- 
come apparent, when the coarsening of domains related 
to the conserved second mode is investigated. 



C. The Gross-Pitaevskii equation 



We finally investigate domain formation in the Gross- 
Pitaevskii equation. As has been remarked earlier, this 
formalism assumes that the dynamics of the condensate 
is completely determined by the precessional (and not re- 
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laxational) dynamics of the classical order parameter. I 
is thus a formalism that on the one hand is valid strictl; 
at zero temperature, but on the other hand ignores quan 
tum fluctuations. The initial state chosen for models I 
and F considered earlier does not evolve in this formal 
ism and hence we choose a slightly different initial stati 
as mentioned in section VII with 90% of the condensati 
density in the state and 5% each in the +1 and -1 states 
The precessional nature of the GP equation implies tha S 
there is never any true relaxation to a state with onh 
domains of +1 and -1 and the amplitudes of these twi 
components oscillate together, 7r/2 out of phase with thi 
amplitude of the component. The dynamical critica 
exponent z in these simulations is extracted by lookinj 
at the time interval when the amplitudes of the +1 ant 
-1 components are growing with time and that of the ( 
component falling. It is important that a sizeable win 
dow be identified within this interval where the domaii 
size L(t) is indeed growing as L{t) oc The oscil 

latory nature of the dynamics ensures that in this case, 
both the magnetization and condensate density are con- 
served. The former is manifested in the fact the +1 and 
-1 components always have the same amplitude and the 
latter in the fact that these two components are always 
7r/2 out of phase with the component. 
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FIG. 9: llzif) as a function of 1/L{t) in the GP equation 
witli the parameter set TZ. The exponent z{t) as a function 
of t and thus increasing L{t) seems to drift towards the value 
z — 3 like with models A and F with conserved magnetization 
densities. 



l°9io* 



FIG. 8: L{t) as a function of log^Q t in the GP equation with 
the parameter set TZ. The data is from a time interval during 
which the amplitudes of the ±1 components are increasing 
with time. 

The domain size here is obtained only using Eqn. 1361 
since the presence of bubbles of the state renders the 
method of measuring the domain boundaries directly un- 
reliable. The data for the domain size L(t) as a function 
of t is shown in Fig. [H The data as in the previous 
two cases has been obtained over a range of about four 
decades. The dynamic critical exponent z{t) as extracted 
is shown as a function of 1/L{t). Once again, there ap- 



pears to be a drift towards the value z = 3 at infinite 
time, although in this case, it appears (to within the 
noise) that the drift is more consistent with Eqn. [35] (i.e. 
linear in 1/L{t)) than for models A and F. 

It appears that the GP model gives a different dynamic 
critical exponent z = 3 from models A and F (z = 2), 
unless magnetization is conserved explicitly in the latter. 
This should be compared and contrasted with the case 
of spinless bosons, where it has been hypothesized that 
model F and the GP equation give the same dynamic 
critical exponent for phase ordering^^. Further, this ex- 
ponent was numerically found to be equal to 1 (different 
from model A, which has z = 2) from numerical studies 
of the GP equation in this case. The situation we ana- 
lyze is different from the case of spinless bosons in that 
we study the magnetization. The second sound mode of 
model F arises from total energy and number conserva- 
tion between the condensate and the "normal state" . The 
GP equation too conserves both these quantities. How- 
ever, for bosons with spin, the GP equation also conserves 
magnetization, which is not present in model F unless in- 
cluded by hand. Thus the value of z for the GP equation 
is different from that of model F and agreement is ob- 
tained only when magnetization is explicitly conserved 
in the latter. It is interesting to note however that the 
value of z obtained from the GP equation is larger than 
from model A in our study whereas for spinless bosons it 
is smaller. 
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IX. CONCLUSIONS 



To conclude, we have studied the statics and dynamics 
of spin-1 condensates at finite temperature in the pres- 
ence of a magnetic field. We have obtained a ground 
state phase diagram for this system and have focussed 
on the phase that is most amenable to a numerical study 
of magnetic domain formation in 2D, namely the Ferro- 
magnetic out-of-plane phase and have numerically deter- 
mined the nature and order of the superfiuid and mag- 
netic phase transitions. We have argued that the "cor- 
rect" dynamical model for spinor condensates at finite 
field and temperature is model F in the Halpering and 
Hohenberg classification and have demonstrated that this 
model contains all the modes in the standard GP equa- 
tion. We have numerically studied magnetic domain for- 
mation in the GP model and models A and F with and 
without magnetization conservation, and have found that 
it is the only when magnetization is explicitly conserved 
in models A and F that the dynamic critical exponenent 
z obtained from the GP equations agrees with the z ob- 
tained from them. In the absence of this conservation, 
models A and F yield z — 2. 

While we have focussed only on one part of the phase 
diagram of spinor condensates in this study to highlight 
the difference between various dynamical models, similar 
studies can be performed in the other parts of the phase 
diagram as well. This will be reported elsewhere. It is our 
belief that model F is fundamentally a more complete dy- 
namical model to describe spinor condensates than the 
GP model. In addition to studies of the dynamics far 
away from the critical point, as presented here, this dy- 
namical model could be used to obtain dynamical criti- 
cal exponents, for comparison to dynamical experiments 
near the static phase transition of the spinor condensate. 



APPENDIX A: RELATION OF THE MODEL F 
PARAMETERS TO MEASURABLE QUANTITIES 



As remarked earlier. Model F has many more param- 
eters than the GP equation. Here we comment on how 
these parameters can be related to experimentally mea- 
surable quantities. There are 9 important real param- 
eters, which are Fi, g^, Aq, /i, cq, C2, C, 7. Some 
of them are directly measurable. For example, Aq is the 
thermal conductivity, and C is the specific heat. We now 
discuss how to measure all the other parameters. 



1- ffo 



go only exists in the dynamical equations and does not 
appear in the static free energy. In the general formalism, 
at the operator level, it appears in the Poisson bracket 
between m and as 

[^l,m]=5oV'l (Al) 

This means if we create a particle through , the expec- 
tation value of m in the system increases by g^. Since m 
is effectively the "heat" in the system (this can be seen 
from the physical meaning of C and A), g^ is effectively 
the "heat" per particle. Writing, 

50 = Ta, (A2) 

where is the entropy per particle, the value of go can 
now be obtained from measuring the specific heat 

go = T£dT^. (A3) 
Here c is the specific heat per particle. 



2. Co + Cj^ 



The reason we discuss cq and 7 together is that cq 
and 7 can only appear in the combination cq + 7^ in 
all static quantities. This can be seen from mean field 
theory: m appears in a quadratic term l/(2C)m^ and a 
linear term 7771710: the expectation value of m is thus C7. 
Plugging this value of m into Eqn. I13i the m'^ interaction 
term becomes cq + Cj^, which we denote as c'. How 
do we measure c'? At low temperature, almost all the 
atoms are in the condensate. All the terms in Eqn.[T51are 
proportional to the density of the condensate, except the 
c' term which is proportional to the square of the density. 
This term will therefore contribute to the compressibility 
K of the condensate, where 

When the temperature is low enough, the dominant con- 
tribution to the compressibility will be from the conden- 
sate. By measuring the compressibility, we can measure 
the parameter c'. 



3. 



Knowing the value of c', the value of p. is quite straight- 
forward to measure. It can be related to the density of 
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atoms within the condensate, using the relation 



no 



(A5) 



oscillation frequency can be shown to be (1 + \/5)r2C2no. 
If we know the value of C2 from static experiments, we 
can obtain the value of 



C2 



6. Co and 7 



C2 is a static parameter and should be measurable 
from static properties, for example, the spin suscepti- 
bility. The terms in the free energy involving the mag- 
netization can be rewritten as C2Af| -|- hMz, where h is 
the magnetic field along z. The spin susceptibility is ap- 
proximately h/c2, from which we can deduce 02- In most 
practical cases at low temperature, the value of C2 is not 
very different from that obtained from the atomic s wave 
scattering rates. 



We have shown that cq + can be determined by 
measuring c'. To disentangle the two quantities, we look 
at the second sound velocity. The equation for the second 
sound velocity is 



C 



5o7^ 



2ri(c' + ^)no 
i 2 



(A6) 



If we know Cg, the only unknown variable is 7. Having 
obtained 7 from this equation, we can calculate cq, from 
the known value of c'. 



The mode </)+i + 4>-i will have an oscillatory compo- 
nent and also a part that is decaying. The value of the 
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